#-----------------------------------------------------------------------------------------------------------------
#------------------------------------------Appendix Tables and Figures--------------------------------------------
#-----------------------------------------------------------------------------------------------------------------

#R version 4.0.3 (2020-10-10)

## Required packages (May need to restart R program after installation)
install.packages("xtable")
install.packages("estimatr")
install.packages("multcomp") # may require "mvtnorm", "survival", "TH.data", "MASS" packages
install.packages("readstata13")
install.packages("plyr")
install.packages("dplyr")
install.packages("stargazer")

## Load libraries
library(xtable) # for producing tables
library(plyr) # for data processing
library(dplyr) # for data processing
library(estimatr) # for estimation
library(multcomp) # for general linear hypotheses testing
library(readstata13) # for loading data
library(stargazer) # for producing tables

rm(list = ls())

## Set working directory
#Please set the working directory to the folder where the readme.txt file is located
setwd("xxx/xxx/xxx")

## load college survey data
survey2nd <-read.dta13("data/collegesurvey.dta")

## generate factor treatment variable
survey2nd$groupc_fct <- factor(survey2nd$qc_group, labels = c ("Control", "Social", "Political", "Social.Political"))

## create a file to standardize variables
survey2nd_std <- survey2nd


## generate a social-order treatment dummy for factorial regression
survey2nd_std <- survey2nd_std %>% 
  mutate(treatcvorg = treat_cv) %>%
  as.data.frame()

survey2nd_std <- survey2nd_std %>%
  mutate(treat_cv = replace(treat_cv, treat_cv_po==1, 1) ) %>%
  as.data.frame()

#table(survey2nd_std$treat_cv)
#table(survey2nd_std$treatcvorg)

## generate a repression treatment dummy for factorial regression
survey2nd_std <- survey2nd_std %>% 
  mutate(treatpoorg = treat_po) %>%
  as.data.frame()

survey2nd_std <- survey2nd_std %>%
  mutate(treat_po = replace(treat_po, treat_cv_po==1, 1) ) %>%
  as.data.frame()

#table(survey2nd_std$treat_po)
#table(survey2nd_std$treatpoorg)

## standardize variables
survey2nd_std$scsapproval <- scale(survey2nd_std$sc_legit)
survey2nd_std$distrust <- scale(survey2nd_std$op_takeother)
survey2nd_std$distrustexp <- scale(survey2nd_std$op_takeself)
survey2nd_std$propaganda <- scale(survey2nd_std$scinfo_state)
survey2nd_std$femalestd <- scale(survey2nd_std$female)
survey2nd_std$agestd <- scale(survey2nd_std$age)
survey2nd_std$incomestd <- scale(survey2nd_std$faminc)
survey2nd_std$ccpstd <- scale(survey2nd_std$party_status)
survey2nd_std$mediastd <- scale(survey2nd_std$media_phone)
survey2nd_std$behaviorstd <- scale(survey2nd_std$behavior)
survey2nd_std$servicestd <- scale(survey2nd_std$par_org)
survey2nd_std$scoreselfstd <- scale(survey2nd_std$sc_selfscore)

survey2nd_std$treatcvstd <- scale(survey2nd_std$treat_cv)
survey2nd_std$treatpostd <- scale(survey2nd_std$treat_po)
survey2nd_std$treatcvpostd <- scale(survey2nd_std$treat_cv_po)


## generate a variable indicating state media as the only source of information about the SCS
survey2nd_std$infostateonly <- survey2nd_std$scinfo_state

survey2nd_std <- survey2nd_std %>%
  mutate(infostateonly = replace(infostateonly, scinfo_nonstate==1, 0) ) %>%
  as.data.frame()

# check the variable
table(survey2nd_std$infostateonly)
#  0   1
#557 180

## load nationwide survey data
nsurvey <- read.dta13("data/nationwidesurvey.dta")

## generate variable female
nsurvey$female = 1- nsurvey$gender

## standardize variables
nsurvey_std <- nsurvey

nsurvey_std$scsapproval <- scale(nsurvey_std$credit_systems_approval)
nsurvey_std$distrust <- scale(nsurvey_std$are_chinese_distrusting)
nsurvey_std$avoidfriend <- scale(nsurvey$lower_friends)
nsurvey_std$infotvnewspaper <- scale(nsurvey$infotvnews)
nsurvey_std$infotvnewscmc <- scale(nsurvey$credit_information_TV)
nsurvey_std$femalestd <- scale(nsurvey$female)
nsurvey_std$agestd <- scale(nsurvey$age1)
nsurvey_std$incomestd <- scale(nsurvey$income)
nsurvey_std$ccpstd <- scale(nsurvey$ccp_memberdm)
nsurvey_std$pubemploystd <- scale(nsurvey$pubemploy)
nsurvey_std$pilotstd <- scale(nsurvey$living_in_pilot_city)
nsurvey_std$influencestd <- scale(nsurvey$credit_decision_influence)
nsurvey_std$gtruststd <- scale(nsurvey$government_confidence)
nsurvey_std$fairstd <- scale(nsurvey$credit_score_fairness)
nsurvey_std$edustd <- scale(nsurvey$education)
nsurvey_std$cityruralstd <- scale(nsurvey$city_rural)
nsurvey_std$privacyusestd <- scale(nsurvey$privacyuse)

## function for word wrap in figures
wordwrap<-function(x,len) paste(strwrap(x,width=len),collapse="\n") 

#----------------Table A.2 Nonresponse Rate by Group----------------
nonresponse <- rbind(as.character(round(table(survey2nd$groupc_fct),0)),
                     as.character(round(with(subset(survey2nd, is.na(sc_legit)), table(groupc_fct) ),0)),
                     as.character(round(with(subset(survey2nd, is.na(sc_legit)), table(groupc_fct) )/table(survey2nd$groupc_fct),3))
)
colnames(nonresponse) <- c("Control", "Social", "Political", "Social & Political")
rownames(nonresponse) <- c("N of Obs.", "N of Nonresponses", "Nonresponse Rate")
print(xtable(nonresponse),
      include.rownames = TRUE, include.colnames = TRUE, file = "outputapp/nonresponse.txt")


#------------------------Table A.3: Summary Statistics (Nationwide Online Survey)-------------------

##Summary Statistics
myvarnsurvey <- c("credit_systems_approval", "are_chinese_distrusting", "lower_friends", "infotvnews", 
                  "credit_information_TV", "female", "age1", "income", "education", "city_rural", "ccp_memberdm", "pubemploy", 
                  "living_in_pilot_city", "credit_decision_influence", "government_confidence", 
                  "credit_score_fairness", "privacyuse")

stargazer(nsurvey[myvarnsurvey], title="Summary Statistics (Nationwide Survey)", digits=2, out = "outputapp/stat1.txt")



#----------------Table A.4 Balance Check (College Field Survey)----------------
surveybal <- subset(survey2nd,complete.cases(qc_group))
nrow(surveybal)
#[1] 747

balancevar <- c("age", "female", "faminc", "faminc_sat", "party_status", "par_political", "par_org", "par_community",
                "par_discuss", "media_news", "media_tv", "media_phone", "op_takeother", "par_comment")

surveybal$qc_group_fct <- factor(surveybal$qc_group, labels = c ("Control", "Social", "Political","Social&Political"))


aovfunc <- function(x){
  p_value <- summary(aov(x ~ qc_group_fct, data = surveybal))[[1]][["Pr(>F)"]][1]
  return(p_value)
}

meanfunc <- function(x){var_mean <- mean(x, na.rm = TRUE)
return(var_mean)
}

sdfunc <- function(x){var_sd <- sd(x, na.rm = TRUE)
return(var_sd)
}

countfunc <- function(x){var_cnt <- length(which(!is.na(x)))
return(var_cnt)
}

out.means <- t(surveybal %>% 
                 group_by(qc_group_fct) %>% 
                 summarise_each (funs(meanfunc),balancevar))

out.pvals <- t(surveybal %>% 
                 summarise_each (funs(aovfunc),balancevar))

out.cnts <- t(surveybal %>% 
                summarise_each (funs(countfunc),balancevar))

out.tab2 <- cbind(out.cnts, apply(out.means[-1,],2, as.numeric),out.pvals)

rownames(out.tab2) = c("Age", "Female (F=1)", "Income (1-9)", "Income Sat. (0-10)", "Party (Yes=1)", "Official Organization (Yes=1)",
                       "Student Organization (Yes=1)", "Community Service (1-5)", "Speech (1-5)", "Media: News (1-5)",
                       "Media: TV (1-5)", "Media: Phone (1-5)", "Distrust (0-10)", "Discuss Politics (1-5)")
colnames(out.tab2) = c("obs.", "Control", "Social", "Political", "Social&Political", "p-value")


print(xtable(out.tab2, digits = c(0,0,2,2,2,2,3)),
      include.rownames = TRUE, include.colnames = TRUE, file = "outputapp/stat2.txt")


#------------------------Figure B.1: Information Treatment Effects, By Online Activism-------------------

survey2nd_std$polcommentonline <- survey2nd_std$par_comment

survey2nd_std <- survey2nd_std %>%
  mutate(polcommentonline = replace(polcommentonline, par_comment==1, 0) ) %>%
  mutate(polcommentonline = replace(polcommentonline, par_comment==2, 0) ) %>%
  mutate(polcommentonline = replace(polcommentonline, par_comment==3, 1) ) %>%
  mutate(polcommentonline = replace(polcommentonline, par_comment==4, 1) ) %>%
  mutate(polcommentonline = replace(polcommentonline, par_comment==5, 1) ) %>%
  as.data.frame()

table(survey2nd_std$polcommentonline)
#  0   1 
#478 263 

mod.expnonsal <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd + factor(univ), 
                           data = subset(survey2nd_std, polcommentonline == 0), se_type = "stata")

mod.expsal <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd + factor(univ), 
                        data = subset(survey2nd_std, polcommentonline == 1), se_type = "stata")


testponsal <- summary(glht(mod.expnonsal, linfct = c("treatpostd+treatcvpostd = 0")))
testcvnsal <- summary(glht(mod.expnonsal, linfct = c("treatcvstd+treatcvpostd = 0")))

testposal <- summary(glht(mod.expsal, linfct = c("treatpostd+treatcvpostd = 0")))
testcvsal <- summary(glht(mod.expsal, linfct = c("treatcvstd+treatcvpostd = 0")))

tab.expnonsal0 <- rbind(
  cbind(testponsal$test$coefficients, testponsal$test$coefficients - 1.96* testponsal$test$sigma, testponsal$test$coefficients + 1.96* testponsal$test$sigma),
  cbind(testcvnsal$test$coefficients, testcvnsal$test$coefficients - 1.96* testcvnsal$test$sigma, testcvnsal$test$coefficients + 1.96* testcvnsal$test$sigma)
)

tab.expsal0 <- rbind(
  cbind(testposal$test$coefficients, testposal$test$coefficients - 1.96* testposal$test$sigma, testposal$test$coefficients + 1.96* testposal$test$sigma),
  cbind(testcvsal$test$coefficients, testcvsal$test$coefficients - 1.96* testcvsal$test$sigma, testcvsal$test$coefficients + 1.96* testcvsal$test$sigma)
)


tab.expnonsal <- cbind(mod.expnonsal$coefficients[4:2],mod.expnonsal$conf.low[4:2],mod.expnonsal$conf.high[4:2])
tab.expsal <- cbind(mod.expsal$coefficients[4:2],mod.expsal$conf.low[4:2],mod.expsal$conf.high[4:2])

nas<- cbind(NA,NA,NA)

tab.expnonsal <-rbind(tab.expnonsal0, nas, tab.expnonsal)
tab.expsal <-rbind(tab.expsal0, nas, tab.expsal)


pdf("outputapp/expttlsalcomp.pdf", width=6.5, height=5) 
par(mar=c(3.9, 9.5, 1, 1) + 1)
plot(tab.expnonsal[,1],1:6,xlim=c(-0.6,0.6),ylim=c(0.41, 6.41),main = "",ylab = "",pch="",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.expnonsal[1:2,1],1:2, pch=1)
points(tab.expnonsal[4:6,1],4:6, pch=16, cex=1.5)
axis(2,at=c(6,5,4,2,1),labels=sapply(c("Social Order", "Political Repression", "Order & Repression", "Order+Interaction", "Repression+Interaction"
),wordwrap, len=20),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
abline(h=3)
legend(-0.655, 6.65, legend="Point Estimates:", box.lty=0)
legend(-0.655, 2.85, legend="Linear Combinaions:", box.lty=0)
for(i in 1:2) segments(tab.expnonsal[i,2],i,tab.expnonsal[i,3],i,lwd=1)
for(i in 4:6) segments(tab.expnonsal[i,2],i,tab.expnonsal[i,3],i,lwd=1.5)
par(new=TRUE)
plot(tab.expsal[,1],1:6,xlim=c(-0.6,0.6),ylim=c(0.59, 6.59),main = "",ylab = "",pch="",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.expsal[1:2,1],1:2, pch=2)
points(tab.expsal[4:6,1],4:6, pch=17, cex=1.5)
for(i in 1:2) segments(tab.expsal[i,2],i,tab.expsal[i,3],i,lwd=1)
for(i in 4:6) segments(tab.expsal[i,2],i,tab.expsal[i,3],i,lwd=1.5)
legend(0.15, 1.15, legend=c("Inactive Online", "Active Online"), pch=1:2, lty=1, lwd=0.75, cex=0.75,pt.cex=0.75, box.lty=0)
legend(0.15, 4.0, legend=c("Inactive Online", "Active Online"), pch=16:17, lty=1, lwd=1, cex=0.75,pt.cex=1, box.lty=0)
dev.off()


#--------------------------Figure B.2: Information Treatment Effects, No Foreigners and Insincere Survey Takers------------------------
## subset data, drop foreigners and insincere survey takers
survey2ndsubset_std <- subset(survey2nd_std, provcppcc < 1 & id != 72 & provcppcc != 1 & born_prov != "TAIWAN" & born_prov != "CAMBODIA" & born_prov != "RUSSIA") # drop respondents who failed the attention check question.

modsubset.approvalexp <- lm_robust(scsapproval ~ distrust + propaganda + distrustexp +
                                     agestd + femalestd + incomestd + ccpstd + behaviorstd +servicestd + factor(univ), cluster = univ,
                                   data = survey2ndsubset_std, se_type = "stata")
tabsubset.approvalexp <- cbind(modsubset.approvalexp$coefficients[10:2],modsubset.approvalexp$conf.low[10:2],modsubset.approvalexp$conf.high[10:2])

modsubset.exp <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd +  factor(univ),
                           data = survey2ndsubset_std, se_type = "stata")
modsubset.exp$nobs
#[1] 721
modsubset.expinfo <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd + factor(univ), 
                               data = subset(survey2ndsubset_std, infostateonly == 0), se_type = "stata")
modsubset.expinfo$nobs
#[1] 545
modsubset.expnoinfo <- lm_robust(scsapproval ~ treatcvstd + treatpostd + treatcvpostd + factor(univ), 
                                 data = subset(survey2ndsubset_std, infostateonly == 1), se_type = "stata")
modsubset.expnoinfo$nobs
#[1] 176

testpo <- summary(glht(modsubset.exp, linfct = c("treatpostd+treatcvpostd = 0")))
testcv <- summary(glht(modsubset.exp, linfct = c("treatcvstd+treatcvpostd = 0")))

testpoi <- summary(glht(modsubset.expinfo, linfct = c("treatpostd+treatcvpostd = 0")))
testcvi <- summary(glht(modsubset.expinfo, linfct = c("treatcvstd+treatcvpostd = 0")))

testponi <- summary(glht(modsubset.expnoinfo, linfct = c("treatpostd+treatcvpostd = 0")))
testcvni <- summary(glht(modsubset.expnoinfo, linfct = c("treatcvstd+treatcvpostd = 0")))

tabsubset.exp0 <- rbind(
  cbind(testpo$test$coefficients, testpo$test$coefficients - 1.96* testpo$test$sigma, testpo$test$coefficients + 1.96* testpo$test$sigma),
  cbind(testcv$test$coefficients, testcv$test$coefficients - 1.96* testcv$test$sigma, testcv$test$coefficients + 1.96* testcv$test$sigma)
)

tabsubset.expinfo0 <- rbind(
  cbind(testpoi$test$coefficients, testpoi$test$coefficients - 1.96* testpoi$test$sigma, testpoi$test$coefficients + 1.96* testpoi$test$sigma),
  cbind(testcvi$test$coefficients, testcvi$test$coefficients - 1.96* testcvi$test$sigma, testcvi$test$coefficients + 1.96* testcvi$test$sigma)
)

tabsubset.expnoinfo0 <- rbind(
  cbind(testponi$test$coefficients, testponi$test$coefficients - 1.96* testponi$test$sigma, testponi$test$coefficients + 1.96* testponi$test$sigma),
  cbind(testcvni$test$coefficients, testcvni$test$coefficients - 1.96* testcvni$test$sigma, testcvni$test$coefficients + 1.96* testcvni$test$sigma)
)

tabsubset.exp <- cbind(modsubset.exp$coefficients[4:2],modsubset.exp$conf.low[4:2],modsubset.exp$conf.high[4:2])
tabsubset.expinfo <- cbind(modsubset.expinfo$coefficients[4:2],modsubset.expinfo$conf.low[4:2],modsubset.expinfo$conf.high[4:2])
tabsubset.expnoinfo <- cbind(modsubset.expnoinfo$coefficients[4:2],modsubset.expnoinfo$conf.low[4:2],modsubset.expnoinfo$conf.high[4:2])


nas<- cbind(NA,NA,NA)

tabsubset.exp <-rbind(tabsubset.exp0, nas, tabsubset.exp)
tabsubset.expinfo <-rbind(tabsubset.expinfo0, nas, tabsubset.expinfo)
tabsubset.expnoinfo <-rbind(tabsubset.expnoinfo0, nas, tabsubset.expnoinfo)


pdf("outputapp/expttlsubset.pdf", width=6.5, height=5) 
par(mar=c(3.9, 9.5, 1, 1) + 1)
plot(tabsubset.exp[,1],1:6,xlim=c(-0.6,0.6),ylim=c(0.5, 6.5),main = "",ylab = "",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tabsubset.exp[1:2,1],1:2, pch=1)
points(tabsubset.exp[4:6,1],4:6, pch=16, cex=1.5)
axis(2,at=c(6,5,4,2,1),labels=sapply(c("Social Order", "Political Repression", "Order & Repression", "Order+Interaction", "Repression+Interaction"
),wordwrap, len=20),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
abline(h=3)
legend(-0.655, 6.65, legend="Point Estimates:", box.lty=0)
legend(-0.655, 2.85, legend="Linear Combinaions:", box.lty=0)
for(i in 1:2) segments(tabsubset.exp[i,2],i,tabsubset.exp[i,3],i,lwd=1)
for(i in 4:6) segments(tabsubset.exp[i,2],i,tabsubset.exp[i,3],i,lwd=1.5)
dev.off()



#--------------------------Figure B.3: Using State Media for Information about the SCS------------------------
mod.statemedia <- lm_robust(infotvnewspaper ~ avoidfriend + distrust + 
                              agestd + edustd + femalestd + incomestd + pubemploystd + pilotstd +
                              cityruralstd + ccpstd + factor(region2), cluster = region,
                            data = nsurvey_std, se_type = "stata")
mod.statemedia$nobs # check sample size
#[1] 1895
tab.statemedia <- cbind(mod.statemedia$coefficients[11:2],mod.statemedia$conf.low[11:2],mod.statemedia$conf.high[11:2])


pdf("outputapp/statemedia.pdf", width=7.9, height=5) 
par(mar=c(3.9, 15, 1, 1) + 1)
plot(tab.statemedia[,1],1:10,xlim=c(-0.335,0.335),ylim=c(0.5, 10.5),main = "",ylab = "",pch=16,cex=1, cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
axis(2,10:1,labels=sapply(c("Avoid Friends with Bad Credits","Perceived Distrust","Age", "Education", "Female", "Income", "State Employee", "Pilot City", "Urban", "CCP Member"
),wordwrap, len=36),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
for(i in 1:9) segments(tab.statemedia[i,2],i,tab.statemedia[i,3],i,lwd=1)
for(i in 9:10) segments(tab.statemedia[i,2],i,tab.statemedia[i,3],i,lwd=1)
dev.off()


#--------------------------Figure B.4: Support for SCSs: Nationwide Online Survey, Commercial Users Only------------------------
mod.approval3 <- lm_robust(scsapproval ~ infotvnewscmc + avoidfriend + distrust + 
                             agestd + edustd + femalestd + incomestd + pubemploystd + pilotstd +
                             cityruralstd + ccpstd + factor(region2), cluster = region,
                           data = nsurvey_std, se_type = "stata")
mod.approval3$nobs # check sample size
#[1] 1410
tab.approval3 <- cbind(mod.approval3$coefficients[12:2],mod.approval3$conf.low[12:2],mod.approval3$conf.high[12:2])


pdf("outputapp/approval3.pdf", width=7.9, height=5) 
par(mar=c(3.9, 15, 1, 1) + 1)
plot(tab.approval3[,1],1:11,xlim=c(-0.335,0.335),ylim=c(0.5, 11.5),main = "",ylab = "",pch="",cex.lab=1,cex.axis=1, yaxt="n",xlab="Coefficients of Standardized Variables with 95% CIs")
points(tab.approval3[1:9,1],1:9, pch=1)
points(tab.approval3[10:11,1],10:11, pch=16, cex=1.5)
axis(2,11:1,labels=sapply(c("SCS Info. Source: TV & Newspaper", "Avoid Friends with Bad Credits","Perceived Distrust","Age", "Education", "Female", "Income", "State Employee", "Pilot City", "Urban", "CCP Member"
),wordwrap, len=36),cex.axis=1,mgp=c(2,1,0),las=1)
#mtext("Treatment",side=1,line=6.5,cex=1.5)
abline(v=0,lty=2,col="red")
for(i in 1:9) segments(tab.approval3[i,2],i,tab.approval3[i,3],i,lwd=1)
for(i in 10:11) segments(tab.approval3[i,2],i,tab.approval3[i,3],i,lwd=1.5)
dev.off()


